NAT'L  INST.  OF  STAND  & TECH  R.I.C. 

Illlll  ' 


AlllDM  DblEat. 


Applied  and 
Computational 
Mathematics 
Division 


Computing  and  Applied  Mathematics  Laboratory 


Boundary/Interface  Fitted  Grid  Generation 
Using  Tensor  Product  B-splines:  A 
Preliminary  Study 

B.V.  Saunders 

August  1993 


Technology  Administration 
U.S.  DEPARTMENT  OF  COMMERCE 
National  Institute  of  Standards  and  Technology 
-QQ Gaithersburg,  MD  20899 

100 

.U56 

#5239 

1995 


Boundary/lnterface  Htted  Grid 
Generation  Using  Tensor  Product 
B-spiines:  A Preiiminary  Study 


B.  V.  Saunders 


U.S.  DEPARTMENT  OF  COMMERCE 
Technology  Administration 
National  Institute  of  Standards 
and  Technology 

Applied  and  Computational  Mathematics  Division 
Computing  and  Applied  Mathematics  Laboratory 
Gaithersburg,  MD  20899 

August  1993 


U.S.  DEPARTMENT  OF  COMMERCE 
Ronald  H.  Brown,  Secretary 

TECHNOLOGY  ADMINISTRATION 

Mary  L Good,  Under  Secretary  for  Technology 


NATIONAL  INSTITUTE  OF  STANDARDS 
AND  TECHNOLOGY 
Arati  Prabhakar,  Director 


■ '■’^’■/v ;- ■ If-* -'-■ ' '^> 'Sit '*•■■*■'  ^-.t.t^K'  - ■'  'S--  ^',.'.'i,uMH 


T,.v 


,•  -^.-lih.  ’""V 

'«**  'i'',.',  'tl',.j!-'.  M?. 

'i  . i\  ■ * 

K /V’V 


tis  -- 


%>  14 


' j 

'1:# 


.■  ~ )l , 


, ),(' 


•■'i  i.U.  ■ - 


. (,..  r - ■ '■  " 

I'  !' 

, ■ - . . n 

'.  ’.  tii,  <i  virii't 

•l-'  -..t,!..  !..^7 

";"■  ■ , ■■  ■ ( <>,11  .„  ;■ 


'*s«.  -■■  : ..  . . 

. . ■ifji 


f-  ,i 

:'  ...y 


% ' # 


•;r^^' 


"ij 


Uj.  V .. 


vT-'V 


Boundary/Interface  Fitted  Grid  Generation  Using 
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National  Institute  of  Standards  and  Technology 
Gaithersburg,  MD  20899 
USA 


Abstract 

Progress  in  the  development  of  an  algebraic  grid  generation  system  that  tracks  a 
solid-liquid  interface  during  directional  solidification  of  a binary  alloy  is  discussed.  A 
single  mapping,  constructed  with  tensor  product  B-splines,  is  proposed  for  calculations 
of  both  shallow  and  deep  solidification  cells.  The  initial  spline  coefficients  for  the 
coordinate  mapping  are  modified  to  minimize  a discrete  functional  that  regulates  the 
smoothness  and  orthogonality  of  the  mesh.  The  use  of  transfinite  blending  function 
interpolation  to  obtain  an  initial  grid  is  examined. 

Keywords:  boundary  fitted  grid  generation,  algebraic  grid  generation,  adaptive,  B- 
splines,  transfinite  blending  functions,  directional  solidification 


1.  Introduction 

One  of  the  well  known  techniques  used  to  study  the  microstructures  that  develop  during 
sohdification  of  binary  alloys  is  Bridgman  growth,  a directional  solidification  technique  in 
which  a sample  of  the  alloy  is  drawn  through  a constant  temperature  gradient  at  a uniform 
rate  of  speed,  V,  as  shown  in  Figure  1.  A considerable  amount  of  theoretical  work  has 
focused  on  examining  the  morphological  instabilities  of  the  growing  soHd-liquid  interface  [1- 
6].  MulHns  and  Sekerka  [1]  used  Hnear  stability  theory  to  predict  the  critical  velocity  for  the 
onset  of  instabihty  for  a planar  interface.  Experimental  observations  confirm  the  vahdity  of 
their  results  and  show  that  after  the  onset  of  instability  the  structure  of  the  interface  can 
change  from  planar  to  cellular  to  dendritic  and  back  again  as  the  growth  velocity  is  increased 
[7,  8].  Coriell  et  al.  [2]  extended  the  results  of  MulHns  and  Sekerka,  including  the  effects 
from  convection  in  the  Hquid. 

Although  cellular  microstructures  are  much  simpler  than  dendritic,  as  the  control  param- 
eters, growth  velocity  or  temperature  gradient,  are  changed,  the  ceUs  may  become  very  deep 
and  narrow  with  re-entrant  bulb-Hke  shapes.  To  successfully  track  the  interface,  the  grid 
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Figure  1:  Bridgman  growth  technique. 

generation  mapping  must  adapt  to  large  deformations  of  the  interface  shape  while  maintain- 
ing as  much  orthogonality  and  smoothness  as  possible.  Ettouney  and  Brown  [4]  successfully 
modeled  slightly  nonplanar  interfaces  by  using  an  algebraic  grid  generation  system  where  the 
interface  was  described  in  terms  of  a univariate  function.  Unfortunately,  this  representation 
of  the  interface,  called  a Monge  transformation,  fails  for  cells  with  vertical  or  re-entrant 
walls.  To  overcome  that  problem,  Ungar  and  Brown  [5]  developed  a representation  defined 
by  the  division  of  the  interface  into  disjoint  sections  that  could  be  expressed  as  separate 
Monge  transformations  written  in  either  polar  or  rectangular  coordinates.  With  this  mixed 
transformation,  Ungar  and  Brown  modeled  cellular  interfaces  with  grooves  as  much  as  15 
times  longer  than  their  wavelength.  Ideally,  analysts  want  a single  grid  generation  tech- 
nique that  can  track  the  interface  as  it  changes  from  a shallow  deformation  to  a narrow  and 
deeply  grooved  cell  with  re-entrant  sidewalls.  Furthermore,  it  should  capture  the  change  in 
wavelength  if  the  cell  sphts. 

Tsiveriotis  and  Brown  [6]  have  made  some  progress  in  this  area  by  using  a two-step 
procedure.  In  one  direction  the  coordinate  is  defined  by  a generaHzed  Poisson  equation  with 
a scahng  condition  that  forces  the  grid  lines  into  concave  areas  of  the  interface.  The  other 
coordinate  is  obtained  by  minimizing  smoothness  and  orthogonality  functionals  similar  to 
those  of  BrackbiU  and  Saltzman  [9].  Tsiveriotis  and  Brown  show  impressive  examples  of 
grids  developed  for  narrow  and  re-entrant  cells,  but  it  is  not  clear  whether  this  method  is 
as  successful  as  the  previous  method  of  Ungar  and  Brown  because  the  most  narrow  cells  in 
the  earher  paper  are  not  redone.  Also,  to  obtain  better  results,  they  fix  a horizontal  line 
near  the  interface  to  decouple  the  far  field  domain  from  the  domain  near  the  interface.  This 
essentially  creates  two  separate  domains  whose  grid  cells  are  not  smoothly  connected  near 
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the  line.  This  is  fine  for  finite  element  calculations,  but  suggests  problems  for  those  who  are 
interested  in  using  the  technique  to  create  grids  suitable  for  finite  difference  schemes. 

This  paper  describes  progress  to  date  in  the  development  of  an  algebraic  technique  for 
generating  a boundary/interface  fitted  coordinate  system  that  is  smooth  enough  to  be  used 
for  finite  difference  calculations.  No  partial  differential  equations  (pdes)  are  solved  to  obtain 
the  coordinate  system,  and  the  same  system  is  used  for  the  entire  domain.  The  computation 
of  the  interface  is  not  discussed  in  this  paper.  It  is  assumed  that  the  interface  is  available  as 
a discrete  set  of  points,  having  been  determined  by  some  means  such  as  an  evaluation  of  the 
Gibbs- Thomson  equation  which  relates  the  interface  mean  curvature,  concentration  of  solute 
in  the  liquid,  and  the  melting  temperature  near  the  interface  [8,  10].  The  grid  generation 
algorithm  is  an  extension  of  an  algebraic  technique  for  generating  boundary  fitted  grids  [11]. 
This  paper  discusses  the  modifications  needed  and  addresses  the  problems  in  tracking  an 
interface  that  deforms  from  a planar  shape  to  a deep,  re-entrant  cellular  microstructure. 

Section  2 provides  a brief  description  of  algebraic  techniques  for  generating  coordinate 
systems,  focusing  on  transfinite  blending  function  interpolation,  the  technique  used  to  obtain 
an  initial  approximation  for  the  boundary /interface  fitted  mesh.  Section  3 explains  the 
construction  of  the  grid  generation  mapping  and  Section  4 examines  the  results  to  date. 


2.  Algebraic  Grid  Generation  Techniques 

In  algebraic  grid  generation,  a direct  transformation  describes  the  relationship  between  the 
computational  and  physical  domains.  The  transformation  is  constructed  so  that  it  interpo- 
lates the  boundary  points  and/or  points  in  the  interior.  It  may  also  be  constrained  to  match 
derivatives.  No  partial  differential  equations  are  solved  to  obtain  the  curvilinear  coordinates, 
so  algebraic  techniques  can  be  easier  to  construct  than  pde  methods,  and  give  easier  con- 
trol over  grid  characteristics  such  as  orthogonality  and  grid  point  spacing.  However,  these 
methods  are  sometimes  criticized  for  allowing  discontinuities  on  the  boundaxy  to  propagate 
into  the  interior  and  for  not  generating  grids  as  smooth  as  those  generated  by  pde  methods. 
Nevertheless,  algebraic  techniques  have  been  used  successfully  to  generate  grids  in  both  two 
and  three  dimensions  [11-16]. 

One  of  the  most  common  and  easiest  algebraic  methods  to  implement  is  transfinite 
blending  function  interpolation,  where  interior  regions  are  represented  in  terms  of  bound- 
ary functions.  This  technique,  which  has  been  widely  used  for  problems  in  grid  generation 
[11,  12,  14,  15,  16],  was  originally  developed  for  problems  in  computer  aided  design.  Gor- 
don and  Hall  [14]  adapted  and  apphed  the  technique  to  grid  generation  for  finite  element 
and  finite  difference  calculations  in  the  late  sixties  and  early  seventies.  A simple  example  is 
illustrated  by  the  mapping  T from  the  unit  square  I2  to  the  physical  domain  defined  by 

T«.r,)  = 

t 
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(1) 


i,3 


where  0 = {o  < • • • < Cm  = 1 and  0 = 7]o  < . . . < r]N  = 1.  The  continuous  vector- 
valued function  f describes  the  boundary  of  the  physical  domain.  $ and  are  “blending” 
or  “connecting”  functions  that  satisfy  the  conditions 

^.(0)  = = I o’  i ^ j (2) 

Mm)  = 4,  = { J;  ^ ; (3) 

Therefore,  if  M=N=1,  one  might  choose  $ and  to  be  hnear  Lagrange  polynomials  so  that 
T matches  f on  the  boundary  of  fy-  If  M and  N are  larger  and  derivative  information  is 
given,  the  blending  functions  can  be  chosen  so  that  interior  curves  and  the  derivatives  there 
are  matched.  Unfortunately,  transfinite  blending  function  interpolation  may  allow  boundary 
singularities  to  propagate  into  the  interior  of  the  mesh.  Furthermore,  if  the  boundary  is 
nonconvex,  the  grid  lines  may  overlap  the  boundary.  The  choice  of  blending  functions  and 
parametrization  of  the  boundary  function  are  very  important  in  alleviating  these  problems. 
Transfinite  blending  function  interpolation  provides  a relatively  easy  way  of  obtaining  an 
initial  grid  that  can  be  refined  and  smoothed  by  other  techniques,  whether  algebraic,  pde, 
or  variational.  However,  unfortunately,  some  techniques  cannot  use  the  grid  if  it  contains 
areas  where  the  grid  lines  overlap,  that  is,  areas  of  negative  Jacobian  [17]. 

Saunders  [11]  developed  a boundary  fitted  grid  generation  mapping  composed  of  tensor 
product  B-splines  whose  coefficients  are  originally  chosen  so  that  the  mapping  approximates 
transfinite  blending  function  interpolation.  An  apphcation  of  the  technique  is  shown  in 
Figures  2 and  3 for  a puzzle  shaped  domain.  The  boundary  and  initial  grid  are  shown 
in  Figure  2.  The  initial  grid  was  constructed  using  hnear  Lagrange  polynomials  for  the 
blending  functions.  Note  that  the  grid  fines  overlap  the  nonconvex  boundary.  The  grid 
generation  system  iterates  on  the  coefficients  to  minimize  a functional  designed  to  increase 
the  smoothness  and  orthogonality  of  the  initial  mesh.  The  grid  on  the  left  in  Figure  3 
shows  the  mesh  obtained  using  the  minimized  coefficients.  The  skewed  areas  have  been 
eliminated  and  the  overlapping  grid  fines  have  been  pulled  into  the  interior.  The  grid  on 
the  right,  computed  with  the  same  coefficients,  shows  a refined  mesh  concentrated  near  the 
bottom  boundary  using  an  exponential  function.  One  dimensional  concentrations  of  this 
type  are  easy  to  obtain  by  replacing  the  mapping  T(fy7/)  with  T(^(^),  7/(77))  where  ( and  fj 
are  appropriate  spacing  control  functions  in  C and  7/  [11].  The  algorithm  discussed  in  this 
paper  is  an  extension  of  the  boundary  fitted  mapping  to  interpolate  a curve  in  the  interior 
of  the  physical  domain.  As  in  the  case  of  the  boimdary  fitted  algorithm,  it  will  be  shown 
that  the  new  algorithm  is  robust  enough  to  smooth  and  untangle  an  initial  grid  containing 
skewed  and  overlapping  grid  lines. 
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Figure  2:  Puzzle  shaped  domain.  Grid  on  right  was  produced  using  an  approximation  of 
transfinite  blending  function  interpolation. 


Figure  3:  Optimized  puzzle  grids. 
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3.  Grid  Generation  Mapping 


The  proposed  grid  generation  mapping,  T,  maps  the  unit  square,  I2,  onto  the  physical 
domain  and  is  constructed  so  that  the  interface  is  the  coordinate  curve  77  = 1 /2  as  shown  in 
the  figure.  The  mapping  is  defined  by 
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Figure  4:  Grid  generation  mapping. 
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where  0 < ^,77  < 1 and  each  Bij  is  the  tensor  product  of  cubic  B-sphnes.  Hence,  = 

Bi{^)Bj{Tj)  where  Bi  and  Bj  are  elements  of  cubic  B-spline  sequences  associated  with  finite 
nondecreasing  knot  sequences,  say,  and  respectively.  As  defined  by  de 

Boor  [18],  B-sphnes  are  essentially  piecewise  polynomials  with  continuity  conditions  at  each 
breakpoint  determined  by  the  repetition  of  the  breakpoint  value  in  the  associated  knot 
sequence.  For  a typical  cubic  B-spline  Bj,  the  value  of  the  B-spHne  is  determined  by  the  five 
knots  tj,  tj+i,  tj+2,  tj+3i  ij+4-  Its  support  is  small,  that  is,  Bj  can  be  nonzero  only  on  the 
interval  [tj,tj+4].  Consequently,  only  four  B-splines,  Bj-z,  Bj_2,  Bj-i,  Bj  can  be  nonzero  on 
the  interval  {tj,tj^i). 

The  spline  coefficients  for  T can  be  divided  into  three  groups.  The  boundary  coefficients 
are  the  coefficients  of  the  Bij  that  are  nonzero  on  the  boundary  of  I2.  Since  T is  defined  so 
that  the  interface  corresponds  to  coordinate  curve  77  = 1/2,  the  coefficients  of  the  Bij  that 
are  nonzero  when  77  = 1/2  are  called  the  interface  coefficients.  The  remaining  coefficients 
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are  called  the  interior  coefficients.  To  obtain  an  initial  grid,  the  sphne  coefficients  are  chosen 
to  produce  a variation  diminishing  spline  approximation  to  the  transfinite  blending  function 
interpolant  that  matches  the  boundary  and  interface  of  the  physical  domain.  In  short,  this 
means  the  coefficients  are  obtained  by  evaluating  the  transfinite  mapping  at  average  knot 
values  as  discussed  in  [18].  This  shape  preserving  approximation  reproduces  straight  lines 
and  preserves  convexity  [11,  18]. 

If  the  deformation  of  the  interface  is  only  moderate,  then  the  initial  grid  may  be  suit- 
able; but  as  the  interface  becomes  more  distorted,  additional  smoothing  and  orthogonahty 
are  needed.  To  increase  orthogonahty  and  decrease  the  overlapping  of  grid  hnes  and  skew- 
ness, the  interior  coefficients  of  T are  modified  to  minimize  a discrete  approximation  to  the 
following  functional  first  described  in  [11]: 


where  T denotes  the  grid  generation  mapping,  J is  the  Jacobian  of  the  mapping,  and  Wi(^,  tj) 
and  W2{C,ri)  are  weight  functions.  When  Wi  is  large,  the  variation  of  the  Jacobian  values  at 
nearby  points  will  be  small,  thereby  decreasing  skewness.  When  W2  is  large,  the  dot  product 
term  will  be  small,  causing  the  grid  lines  to  approach  orthogonahty. 

Clearly,  the  flexibility  of  the  boundary  and  the  interface  coefficients  is  limited.  Using  the 
continuity  properties  of  B-splines  with  repeated  knots  [18],  one  can  show  that  if  the  first 
four  knots  of  and  are  0,  the  last  four  knots  are  1,  and  the  rest  located  in  the 

interval  (0, 1),  the  only  boundary  coefficients  are  aij,^ij  and  for  j = 1, . . . ,n  and 

0^1  j Ai  3-nd  Oin,  Pin  for  i = 1, . . . , m.  Furthermore,  unless  all  the  “interior”  knots  are  clustered 
near  the  center  of  (0, 1),  the  boundary  coefficients  will  have  a minimal  effect  on  the  interior  of 
the  square.  For  example,  aij,  will  only  affect  points  on  the  support  of  Bij,  that  is,  the  narrow 
band  [0,55]  x [tj,tj+4]  illustrated  in  Figure  5,  where  Ss  is  the  first  interior  knot  in 
In  fact  the  area  of  influence  of  ah  the  boundary  coefficients  wiU  be  a narrow  region  along 
the  boundary  of  the  square.  Consequently,  even  if  the  boundary  coefficients  remain  fixed, 
the  smoothing  functional  should  be  able  to  produce  a significant  amount  of  orthogonahty 
and  smoothness  in  the  grid.  Ideally,  of  course,  the  boundary  coefficients  should  also  be 
adjusted  by  the  smoothing  functional,  but  the  boundary  fittedness  requirement  severely 
restricts  their  movement.  For  a general  boundary  of  arbitrary  shape,  small  changes  in  the 
boundary  coefficients  can  destroy  the  shape  preserving  properties  of  the  spline  mapping.  In 
such  a case,  fixing  the  boundary  coefficients  is  probably  the  simplest  solution.  However,  if  the 
boundary  is  rectangular,  as  it  is  for  Bridgman  growth,  more  flexibility  is  possible.  Using  the 
smoothing  functional,  the  boundary  coefficients  can  be  modified  so  that  a reparametrization 
of  the  boundary  mapping  results. 

Since  the  interface  may  be  quite  complex,  very  httle  change  is  permitted  in  the  interface 
coefficients  when  the  other  coefficients  are  adjusted  by  the  minimization  of  the  smoothing 
functional.  Unfortunately,  unlike  the  boundary  coefficients,  the  interface  coefficients  affect 
a significant  number  of  points  on  the  interior  besides  those  mapped  onto  the  interface.  To 
determine  the  interface  coefficients  one  first  finds  I such  that  1/2  G Then  the 
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Figure  5:  Support  of  tensor  product  B-spline  Bij. 

interface  coefficients  will  be  Pik  where  1 < i < m and  I — 3 < k < 1.  These  are  the 
coefficients  of  the  tensor  product  B-splines  that  might  be  nonzero  when  77  = 1/2.  However, 
fixing  these  coefficients  affects  the  mapping  T not  only  on  [0,1]  x [ii,i/+i],  but  also  on 
the  much  larger  band,  [0,1]  x [t/_3,t/+4]  seen  in  Figure  6,  that  is,  the  total  support  of 
the  associated  tensor  product  B-splines.  This  band  can  be  narrowed  significantly  by  using 
a larger  number  of  knots  for  the  t sequence  and  concentrating  some  near  1/2.  However, 
concentrating  the  knots  too  closely  will  affect  the  smoothness  of  the  grid  lines.  On  the 
other  hand,  since  a continuous  B-sphne  is  close  to  zero  near  the  boundary  of  its  support, 
all  interface  coefficients  do  not  have  an  equal  effect  on  the  interface  mapping.  In  particular, 
allowing  the  coefficients  ai^i-3  and  for  1 < i < m to  move  appears  to  have  very  little 

effect  on  the  accuracy  of  the  interface  mapping.  Nevertheless,  at  the  outset,  the  interface 
coefficients  should  be  chosen  to  produce  as  much  orthogonahty  and  smoothness  as  possible. 
As  expected,  the  results  in  the  next  section  indicate  that  the  type  of  blending  functions  used 
in  defining  the  initial  coefficients  is  very  important. 


4.  Results 

For  aU  the  examples  shown,  constants  of  values  1 and  10  are  used  for  the  Jacobian  and  orthog- 
onahty weights,  respectively,  in  the  smoothing  functional.  As  in  [11],  a simple  minimization 
technique,  the  cyclic  coordinate  method,  is  used  to  minimize  a discrete  approximation  to 
the  functional.  The  number  of  spline  coefficients  used  to  define  the  grid  generation  mapping 
is  2 X X Nr,  where  is  the  number  of  B-sphnes  associated  with  sequence  {5},  or  the 
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Figure  6:  Support  of  tensor  product  B-splines  nonzero  when  rj  = 1J2. 


sequence,  and  Nj,  is  the  number  associated  with  sequence  {t},  the  “t/”  sequence.  A 
concentration  of  grid  points  near  the  interface  (coordinate  curve  rj  — 1/2)  is  obtained  by 
replacing  r}  with  the  hyperbolic  sine  function 

, sinh(2c77  — c)  — sinh(— c) 

^ ^ sinh(c)  — sinh(— c) 

where  c is  a constant  that  determines  the  degree  of  concentration.  To  increase  the  accuracy 
of  the  spline  mapping  and  decrease  the  area  affected  by  the  interface  coefficients,  illustrated 
in  Figure  6,  the  sinh  function  is  also  used  to  concentrate  the  rj  knots  near  the  interface. 
This  means  that  the  sequence  {^(i)}  is  used  instead  of  {t}  in  the  definition  of  the  mapping. 
The  values  of  c associated  with  these  two  types  of  concentration  will  be  called  the  mesh 
concentration  and  the  knot  concentration,  respectively.  It  should  also  be  noted  that  the 
term  “initial  grid”  is  used  to  refer  to  any  grid  computed  using  the  initial  spline  coefficients. 

The  first  example  illustrates  the  difficulty  in  extending  the  algorithm  in  [11]  to  develop  a 
grid  that  interpolates  a solid-liquid  interface.  The  bottom  of  the  puzzle  boundary  displayed 
in  Figure  2 is  now  shown  as  an  interface.  This  is  similar  to  the  re-entrant  shape  that 
commonly  appears  in  cellular  microstructures.  = 19  and  Nr,  = 20  for  this  example.  To 
obtain  the  transfinite  approximation,  the  physical  domain  is  divided  into  top  and  bottom 
halves  separated  by  the  interface.  In  each  section  linear  Lagrange  polynomials  are  used  as 
blending  functions  to  create  a mapping  that  interpolates  the  boundary  and  interface.  The 
initial  coefficients  are  chosen  to  create  a variation  diminishing  spline  approximation  to  the 
transfinite  mapping.  A knot  concentration  of  3 is  used  to  define  the  spline  approximation. 
In  Figure  7 the  mesh  on  the  left  shows  a grid  computed  with  the  initial  coefficients.  As  with 
the  original  puzzle  grid  shown  in  Figure  2,  one  can  see  that  the  grid  cells  near  the  re-entrant 
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curve  are  skewed  and  overlap  each  other.  Attempts  to  improve  the  grid  by  minimizing 
the  smoothing  functional  were  unsuccessful.  The  grid  on  the  right  was  computed  using 
coefficients  obtained  after  twenty  iterations  of  the  cyclic  coordinate  minimization  routine  on 
a 35  X 60  initial  grid  with  mesh  concentration  equal  3. 


Figure  7:  Grids  produced  using  biHnear  blending  functions.  Initial  grid  is  on  the  left.  Grid 
on  right  shows  smoothing  routine  unable  to  improve  grid. 

The  use  of  higher  degree  polynomials  for  the  blending  functions  permits  the  creation  of 
mappings  that  not  only  interpolate  the  boundary  and  interface  but  also  match  derivatives 
there.  For  example,  using  cubic  blending  functions  in  the  rj  coordinate,  one  can  interpolate 
the  top  and  bottom  boundaries,  the  interface,  and  the  normal  vectors  at  the  interface. 
However,  even  when  orthogonality  is  obtained  at  the  interface,  polynomials  of  high  degree 
may  produce  a significant  amount  of  skewness  and  overlap  of  grid  cells  elsewhere  in  the  mesh. 
Among  the  polynomials  tested:  Hnear,  quadratic,  cubic,  and  hermite  quintic  (fifth  degree), 
the  hermite  quintic  polynomials  produced  the  least  amount  of  skewness  and  overlap.  For  aU 
the  grids  that  foUow,  hermite  quintic  blending  functions  were  used  for  the  77  coordinate  and 
linear  Lagrange  polynomials  for  the  ^ coordinate.  Unlike  the  previous  example,  the  same 
transfinite  mapping  is  used  over  the  entire  domain.  Boundary  data  input  by  the  user  is  used 
to  define  a continuous  vector  valued  function  f that  maps  the  boundary  of  the  square  to 
the  boundary  of  the  physical  domain  and  the  line  77  = 1/2  to  the  interface.  The  transfinite 
mapping  T is  constructed  to  interpolate  f on  the  boundary  and  match  f and  df/dr}  at 
77  = 1/2.  It  has  the  form 


i=0 
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7=0 

t=0  7=0 

^ di. 

- («) 

where  ^o  = 0,  = 1 and  770  = 0,  771  = 1/2,  772  = 1.  The  hnear  blending  functions  satisfy 

= ^^3  for  i,j  = 0,1 

while  the  hermite  quintic  blending  functions  and  /3  satisfy 

"^riVs)  = ^TS 

KM  = 0 

PiVs)  = 0 

P'iVs)  = Su  for  r,  s = 0, 1, 2 . 

To  force  orthogonality  at  the  interface,  the  components  of  the  derivative 


are  chosen  to  be 


dr} 

drj 


where 


L = 


'aAV  /%■ 

\{d^)  +Ui. 


and  K is  a.  user  defined  orthogonality  constant  that  regulates  the  magnitude  of  the  normal 
vectors  at  the  interface.  The  variation  diminishing  spline  approximation  to  T produces  the 
mesh  shown  on  the  left  in  Figure  8.  As  in  the  first  example,  the  knot  concentration  equals 
3,  = 19  and  Nr,  = 20.  The  orthogonality  constant  A"  is  3.  Although  the  grid  cells 

are  skewed  in  some  areas,  the  grid  is  more  orthogonal  near  the  interface  than  the  previous 
example.  Hence,  the  initial  grid  already  looks  fairly  good  there.  This  is  important  since 
most  of  the  interface  coefficients  remain  fixed  throughout  the  minimization  process.  Using  a 
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Figure  8:  Grids  produced  using  linear  and  hermite  quintic  blending  functions.  Initial  grid  is 
on  the  left.  Grid  on  right  was  produced  after  twenty  iterations  of  smoothing  routine. 

35  X 60  initial  grid  with  mesh  concentration  3,  the  optimization  routine  easily  improves  the 
smoothness  and  general  orthogonality  of  the  grid  cells.  This  is  illustrated  in  the  grid  on  the 
right,  computed  after  twenty  iterations. 

To  illustrate  the  flexibility  of  the  grid  generation  routine,  the  next  six  grids  show  a typical 
deformation  of  an  interface  from  a planar  shape  to  a deep  cell.  The  grids  were  constructed 
independently  of  each  other.  That  is,  each  grid  was  created  by  first  starting  with  boundary 
data  and  interface  data  for  the  specific  interface  shape.  An  initial  grid  was  formed  from 
the  spline  approximation  to  the  transfinite  mapping  using  a knot  concentration  of  3.  For 
all  the  grids,  K = 2^  and  A/,,  = 24.  The  2 x x Nr,  spline  coefficients  were  then 

iterated  on  to  minimize  the  smoothing  functional.  Except  for  the  planar  grid  which  required 
no  iterations,  the  grids  shown  were  produced  after  twenty  iterations  using  a 30  x 60  grid 
with  mesh  concentration  4.  Although  aU  the  grids  are  relatively  smooth,  one  can  see  that 
smoothness  and  orthogonahty  are  harder  to  maintain  as  the  cell  deepens. 

The  final  example  shows  a grid  constructed  for  a deep  re-entrant  cell.  Figures  12  and 
13  show  the  mesh  after  0,  1,  5,  and  10  iterations,  respectively.  The  orthogonahty  constant 
K = 2,  = 18,  Nr,  = 24  and  the  knot  concentration  is  3.  The  minimization  was  done  on  a 

60  X 50  grid  ivith  mesh  concentration  4.  ]^lote  that  the  initial  grid  overlaps  itself  several  tim.es 
in  the  area  underneath  the  bulb-hke  bottom.  The  smoothing  routine  is  able  to  untangle  this 
area  after  only  five  iterations. 
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Figure  9:  Deformation  of  a solid-liquid  interface.  Stages  0 and  1. 
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Figure  11:  Deformation  of  a solid-liquid  interface.  Stages  4 and  5. 


Figure  12:  Re-entrant  interface  shape.  Initial  grid  on  left  shows  skewed  and  overlapping  grid 
cells.  Grid  on  right  produced  after  1 iteration  of  smoothing  routine. 
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Figure  13:  Re-entrant  interface  shape.  Grid  on  left  produced  after  5 iterations  of  smoothing 
routine.  Grid  on  right  produced  after  10  iterations. 

5.  Conclusions 

Progress  in  the  development  of  an  algebraic  grid  generation  system  that  tracks  a sohd-liquid 
interface  has  been  discussed.  The  proposed  grid  generation  mapping  effectively  interpolates 
interior  curves  of  varying  degrees  of  complexity  from  planar  to  deep  re-entrant  grooved 
shapes.  Transfinite  blending  function  interpolation  appears  to  be  a good  method  for  creating 
an  initial  grid  if  the  blending  functions  are  chosen  appropriately.  Additional  study  wiU  be 
continued  in  this  area. 

Currently,  the  grids  are  developed  by  starting  with  the  boundary  and  interface  data  for 
a particular  shape.  Hence,  the  grid  generation  problem  becomes  quite  challenging  for  deeply 
grooved  and  narrow  re-entrant  shapes.  In  an  actual  application,  the  initial  grid  will  be  the 
grid  used  in  the  previous  iteration.  Therefore,  the  interface  change  should  be  more  gradual. 
Hopefully,  starting  with  a previously  developed  grid  wiU  allow  the  system  to  generate  meshes 
for  interfaces  that  are  more  narrow  and  grooved.  This  is  the  current  area  of  study. 

The  next  phase  wiU  be  to  couple  the  grid  generation  algorithm  with  equations  that 
determine  the  interface  shape  so  that  the  system  can  be  used  in  the  numerical  analysis  of 
microstructures  that  develop  during  directional  solidification.  Although  the  system  will  be 
designed  to  track  the  deformation  of  a planar  interface  into  a deep  cell,  it  may  also  be  used 
to  improve  calculations  in  phase  field  models  where  the  interface  is  not  viewed  as  a curve  or 
surface  with  zero  thickness.  Accuracy  can  be  improved  by  concentrating  the  grid  points  in 
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the  area  near  the  interface  even  if  the  interface  is  not  tracked  exactly. 
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